function y= rootlocus_BD(x,p)
    sigma=x.^p;
    switch p
    case 1
        rho=x-1.0;
    case 2
        rho=x.^2-4/3*x+1/3;
        sigma=2/3*sigma;
    case 3
        rho=x.^3-18/11*x.^2+9/11*x-2/11;
        sigma=6/11*sigma;
    case 4
       rho=x.^4-48/25.*x.^3+36/25.*x.^2-16/25.*x+3/25;
       sigma=12/25*sigma;
    end
    y=1.0*rho./sigma;
end
